clear all;
close all;
delete(gcp('nocreate'));
%Define initial step for estimation algorithm;
global usual_delta;
usual_delta = 0.05;
% Current folder for Matlab codes
fp.matlab = pwd;

%%%%%%%%%%;
%Folders;
%%%%%%%%%%;
% Current folder for Matlab codes
common_path = extractBefore(fp.matlab,"\matlab_dir");
fp.paper = sprintf('%s\\paper', common_path);
fp.boot_dir = sprintf('%s\\matlab_dir\\temp_boot', common_path);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%fixed parameters;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%We estimate the model at lambda=0.95;
fp.lambda = 0.95;

%model primitives;
fp.ages =9:12;
fp.periods=length(fp.ages);
fp.percentiles_kid = linspace(5,95,10);
fp.hc_points_kid = length(fp.percentiles_kid);

fp.Min_Time = 0.01;
fp.Max_Time= 1;
fp.inv_points=20;
fp.percentiles_inv= linspace(5,95,fp.inv_points);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%indexes for simulated data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

fp.ind_data.child_id = 1;
fp.ind_data.child_age = 2;
fp.ind_data.n_sim = 3;
fp.ind_data.child_C = 4;
fp.ind_data.child_C_tp1 = 5;
fp.ind_data.inv = 6;
fp.ind_data.peers = 7;
fp.ind_data.peers_log = 8;
fp.ind_data.ParStyle = 9;
fp.ind_data.peers_tp1 = 10;
fp.ind_data.peers_log_tp1 = 11;
fp.ind_data.neighborhood = 12;
fp.ind_data.moved = 13;

fp.ind_data_network.child_id=1;
fp.ind_data_network.child_age=2;
fp.ind_data_network.n_sim=3;
fp.ind_data_network.child_C=4;
fp.ind_data_network.n_friends=5;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Loading Estimated Moments (and keep track of orders of moments with indexes);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

temp_mean_PS_between_age=importdata('mean_PS_between_class.txt');    
fp.mom_data.mean_PS_between_age = temp_mean_PS_between_age;

moments_temp=importdata('moments.txt');    

variance_moments_temp = importdata('variance_moments.txt');    
fp.variance_moments = variance_moments_temp;

%Store moments and define indexes;
fp.mom_data.mean_PS = moments_temp(1);

fp.ind_mom.coef_PS.child_C = 1;
fp.ind_mom.coef_PS.peers   = 2;
fp.mom_data.coef_PS     = moments_temp(2:3);

fp.ind_mom.coef_skills_tp1.child_C    = 1;
fp.ind_mom.coef_skills_tp1.peers      = 2;
fp.ind_mom.coef_skills_tp1.PS         = 3;

fp.mom_data.coef_skills_tp1 = moments_temp(4:6); %Table 11 Column (1) 
fp.mom_data.coef_skills_tp1_ps0 = moments_temp(7:8);  
fp.mom_data.coef_skills_tp1_ps1 = moments_temp(9:10);

fp.ind_mom.coef_peers_tp1.child_C = 1;
fp.ind_mom.coef_peers_tp1.peers   = 2;
fp.ind_mom.coef_peers_tp1.PS      = 3;

fp.mom_data.coef_peers_tp1  = moments_temp(11:13); %Table 2 Column (2) 
fp.mom_data.coef_peers_tp1_ps0  = moments_temp(14:15); 
fp.mom_data.coef_peers_tp1_ps1  = moments_temp(16:17); 

fp.ind_mom.coef_inv.child_C = 1;
fp.ind_mom.coef_inv.peers   = 2;
fp.ind_mom.coef_inv.const   = 3;

fp.mom_data.coef_inv_PS0  = moments_temp(18:19); 
fp.mom_data.coef_inv_PS1  = moments_temp(20:21); 

fp.mom_data.mean_skills = moments_temp(22:25);
fp.mom_data.mean_inv    = moments_temp(26:27);
fp.mom_data.mean_nfriends = moments_temp(28);


fp.ind_mom.coef_PS_longit.child_C = 1;
fp.ind_mom.coef_PS_longit.peers   = 2;

fp.mom_data.PS_longit = moments_temp(29:30);
fp.mom_data.Inv_longit = moments_temp(31:32);

fp.mom_data.mean_PS_between = moments_temp(33:36);

%Store vector of Moments;
fp.data_mom = moments_temp;

%Loading Initial Conditions (Table 4);
initial_conditions = importdata('initial_conditions.txt');    
fp.schools_distrib = initial_conditions;

%Loading Income information of neighborhoods (it affects preferences for parenting in the model);
mean_income_by_neighborhoods = importdata('family_income_neighborhoods.txt');
mean_income_by_neighborhoods = mean_income_by_neighborhoods./10000;
fp.mean_income_by_neighborhoods = mean_income_by_neighborhoods;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Draws of random mumber (by neighborhoods)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

rng(10); %fix seed;
fp.n_rip_sim = 50;                          % # of simulation of the economy;
fp.tot_draws_mcintegration_network = 50 ;   % # of simulation of the network;

%Initial skills;
fp.init_draws{1} = norminv(rand(fp.schools_distrib(1,3),1,fp.n_rip_sim));
fp.init_draws{2} = norminv(rand(fp.schools_distrib(2,3),1,fp.n_rip_sim));
fp.init_draws{3} = norminv(rand(fp.schools_distrib(3,3),1,fp.n_rip_sim));
fp.init_draws{4} = norminv(rand(fp.schools_distrib(4,3),1,fp.n_rip_sim));

%Peer shocks (backward induction)
fp.peers_backward{1} = rand(length(fp.percentiles_kid)*length(fp.percentiles_kid)*fp.inv_points,fp.schools_distrib(1,3),fp.tot_draws_mcintegration_network, fp.periods-1);
fp.peers_backward{2} = rand(length(fp.percentiles_kid)*length(fp.percentiles_kid)*fp.inv_points,fp.schools_distrib(2,3),fp.tot_draws_mcintegration_network, fp.periods-1);
fp.peers_backward{3} = rand(length(fp.percentiles_kid)*length(fp.percentiles_kid)*fp.inv_points,fp.schools_distrib(3,3),fp.tot_draws_mcintegration_network, fp.periods-1);
fp.peers_backward{4} = rand(length(fp.percentiles_kid)*length(fp.percentiles_kid)*fp.inv_points,fp.schools_distrib(4,3),fp.tot_draws_mcintegration_network, fp.periods-1);

%Peer shocks (forward induction)
fp.peers_forward{1} = rand(fp.schools_distrib(1,3),fp.schools_distrib(1,3),fp.n_rip_sim,fp.periods-1);
fp.peers_forward{2} = rand(fp.schools_distrib(2,3),fp.schools_distrib(2,3),fp.n_rip_sim,fp.periods-1);
fp.peers_forward{3} = rand(fp.schools_distrib(3,3),fp.schools_distrib(3,3),fp.n_rip_sim,fp.periods-1);
fp.peers_forward{4} = rand(fp.schools_distrib(4,3),fp.schools_distrib(4,3),fp.n_rip_sim,fp.periods-1);

%Parenting Style shocks
fp.PS_forward{1} = rand(fp.schools_distrib(1,3),fp.n_rip_sim,fp.periods-1);
fp.PS_forward{2} = rand(fp.schools_distrib(2,3),fp.n_rip_sim,fp.periods-1);
fp.PS_forward{3} = rand(fp.schools_distrib(3,3),fp.n_rip_sim,fp.periods-1);
fp.PS_forward{4} = rand(fp.schools_distrib(4,3),fp.n_rip_sim,fp.periods-1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Estimate the model;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

fp.do_estimation   = 0;
fp.do_counter_type = 0; 

if fp.do_estimation == 1

fp.do_counter = 0;            
                                                                                                                        
    load('param0')              
                  
    options_est = optimset;    
    options_est = optimset(options_est,'MaxFunEvals',3000);
    options_est = optimset(options_est,'MaxIter',3000);
    options_est = optimset(options_est,'TolX',1e+0);
    options_est = optimset(options_est,'TolFun',1e-2);
    options_est = optimset(options_est,'UseParallel',1);
    options_est = optimset(options_est,'Display','iter');
    param_est = fminsearch_function(@obj_function, param0 , options_est,fp );
    save('param_est','param_est');
    
    
end

    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Standard Errors;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

do_se  = 0;
if do_se == 1
fp.do_estimation =1;    
fp.do_counter = 0;                

cd(fp.boot_dir)
%load bootstrap distribution of moments;
moments_bootstrapped=importdata('moments_boot.txt');    
cd(fp.matlab)
load('param0')              
                  
param_est_boot = NaN( 100, length(param0));

for b=1:100        
    fp_parfor_boot = fp;
    %Select the "b" realization of the bootstrap distribution of moments;
    fp_parfor_boot.data_mom = moments_bootstrapped(b,:)';

    cd(fp.boot_dir)
    %load bootstrap distribution of initial conditions;
    initial_conditions_boot          = importdata(sprintf('initial_conditions_boot%d.txt',b));    
    mean_income_by_neighborhoods     = importdata(sprintf('family_income_neighborhoods_boot%d.txt',b));    
    cd(fp.matlab)

    %include n of children;
    fp_parfor_boot.schools_distrib = [ initial_conditions_boot'; fp.schools_distrib(:,3)' ]';

    mean_income_by_neighborhoods = mean_income_by_neighborhoods./10000;
    fp_parfor_boot.mean_income_by_neighborhoods = mean_income_by_neighborhoods;
    
    options_est = optimset;    
    options_est = optimset(options_est,'MaxFunEvals',3000);
    options_est = optimset(options_est,'MaxIter',3000);
    options_est = optimset(options_est,'TolX',1e+0);
    options_est = optimset(options_est,'TolFun',1e-2);
    options_est = optimset(options_est,'Display','iter');
    param_est_boot(b,:)=fminsearch_function(@obj_function, param0 , options_est, fp_parfor_boot );

end
  
save('param_est_boot','param_est_boot')

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%POST-ESTIMATION:

%Table 5: Estimated Parameters of the Skill Formation Technology;
%Table 6: Estimated Parent;s Preference Parameters;
%Table 7: Estimated Child;s Preference Parameters;
%Appendix Table D-1: Sample Fit of the Model: Parenting Style;
%Appendix Table D-2: Sample Fit of the Model: Skill Accumulation;
%Appendix Table D-3: Sample Fit of the Model: Peer Skills;
%Appendix Table D-4: Sample Fit of the Model: Parental Investments;
%Appendix Table D-5: Sample Fit: Longitudinal Analysis of Parenting;

%Appendix Figure D-1: Comparative Statics of Friendship Formation;
%Appendix Figure D-2: Perturbation of Model s Parameters and Equilibrium Moments (Panels A-E);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

do_post_estim = 1;
if do_post_estim==1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Simulate Baseline Economy (and print sample-fit in Matlab command window);       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

load('param_est');   
load('param_est_boot')
fp.param_est_boot = param_est_boot;

format short g
fp.do_estimation = 0;
fp.do_counter    = 0;            
obj_eval = obj_function( param_est, fp );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Tables of Estimates:
%Table 5: Estimated Parameters of the Skill Formation Technology
%Table 6: Estimated Parent’s Preference Parameters
%Table 7: Estimated Child’s Preference Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
cd(fp.matlab) 
Estimates_tables(obj_eval.structural_params,fp); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Tables of Sample Fit:
%Appendix Table D-1: Sample Fit of the Model: Parenting Style
%Appendix Table D-2: Sample Fit of the Model: Skill Accumulation
%Appendix Table D-3: Sample Fit of the Model: Peer Skills
%Appendix Table D-4: Sample Fit of the Model: Parental Investments
%Appendix Table D-5: Sample Fit: Longitudinal Analysis of Parenting
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
cd(fp.matlab)
Sample_fit_tables(obj_eval.sim_moments,fp.data_mom,fp);    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Comparative Statics Exercises;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Appendix Figure D-1: Comparative Statics of Friendship Formation;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
cd(fp.matlab)
graphs_friendships_formation(obj_eval.simul_data,param_est,fp);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Appendix Figure D-2: Perturbation of Model s Parameters and Equilibrium Moments (Panels A-E)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
cd(fp.matlab)
graphs_perturb_parameters(param_est,fp)

end %end if do_post_estim = 1
